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ABSTRACT 

A stochastic adhesion model is introduced, with the purpose of describing 
the formation and evolution of mildly nonlinear structures, such as sheets and 
filaments, in the intergalactic medium (IGM), after hydrogen reionization. The 
model is based on replacing the overall force acting on the baryon fluid - as it 
results from the composition of local gravity, pressure gradients and Hubble 
drag - by a mock external force, self-consistently calculated from first-order 
perturbation theory. A small kinematic viscosity term prevents shell-crossing 
on small scales (which arises because of the approximate treatment of pressure 
gradients). The emerging scheme is an extension of the well-known adhesion 
approximation for the dark matter dynamics, from which it only differs by 
the presence of a small-scale 'random' force, characterizing the IGM. Our 
algorithm is the ideal tool to obtain the skeleton of the IGM distribution, 
which is responsible for the structure observed in the low-column density Lya 
forest in the absorption spectra of distant quasars. 
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1 INTRODUCTION 

The analysis of the spectra of distant QSOs blueward of the Lya emission, the so-called Lya forest, has revealed the 
presence of large coherent structures in the cosmic distribution of the neutral hydrogen, left over by the reionization 
process (e.g. Rauch 1998 and references therein). The low-column density (A^ni < 10^*'^ cm~^) Lya forest is thought as 
being generated by unshocked gas in voids and mildly overdense fluctuations of the photoionized intergalactic medium 
(IGM). The IGM is, in turn, expected to trace the underlying dark matter (DM) distribution on large scales, while the 
thermal pressure in the gas smears small-scale DM nonlinearities. According to this picture, the large coherent spatial 
structures in the IGM, which are probed by quasar absorption spectra, reflect the underlying network of filaments and 
sheets in the large-scale DM distribution (e.g. Efstathiou, Schaye & Theuns 2000), the so-called 'cosmic web' (Bond, 
Kofman & Pogosyan 1996). Thanks to this simple picture, a number of semi-analytical, or pseudo-hydrodynamical 
particle-mesh codes have been developed, which all aim at simulating the IGM distribution, with enough resolution 
to predict the main features of the Lya forest (e.g. Bi & Davidsen 1997; Croft et al. 1998; Gnedin & Hui 1996, 1998; 
Hui, Gnedin & Zhang 1997; Petitjean, Miicket & Kates 1995; Gaztahaga & Croft 1999; Meiksin & White 2000; Viel 
et al. 2001). This physical interpretation of the Lya forest is, moreover, largely corroborated by numerical simulations 
based on fully hydrodynamical codes (e.g. Cen et al. 1994; Zhang, Anninos & Norman 1995, 1997; Miralda-Escude 
et al. 1996; Hernquist et al. 1996; Theuns et al. 1998). 

Although the above scenario captures the main physical processes which determine the coarse-grained dynamics of 
the baryons, it neglects a number of fine-grained details which distinguish the baryon distribution from the underlying 
DM one. There has been, recently, growing observational and theoretical interest on the low-column density Lya 
forest, as an important indicator of the state of the Universe in the redshift range 2 — 4. We then believe it has 
become necessary to develop new and more sophisticated techniques, enabling to understand the IGM dynamics and 
accurately simulate its clustering, without resorting to the full machinery of hydro-simulations. The model presented 
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here aims at providing an accurate scheme for the growth and evolution of mildly nonlinear structures in the baryon 
gas, accounting for the detailed thermal history of the IGM. 

The stochastic adhesion model introduced here is based on applying the forced Burgers equation of nonlinear 
diffusion (Forster, Nelson & Stephen 1977; Kardar, Parisi & Zhang 1986; Barabasi & Stanley 1995; E et al. 1997; 
Frisch & Bee 2000 and references therein) to the cosmological framework, where it can be given the general form 

£}^^^u\/M^,r)-\/v{^,r). (1) 

Here x are comoving coordinates, the time variable r is chosen to coincide with the cosmic scale-factor a (although, 
for the pure DM evolution, a better choice would be the linear growth factor of density fluctuations), u = dx/dr is the 
(suitably rescaled) peculiar velocity field, which is here assumed to be irrotational, u — V$, and D/Dt = d/dr + u-V 
denotes the convective time-derivative. The coefficient of kinematical viscosity, u, is always assumed to be small. In 
standard applications of the forced Burgers equation, ry is a random external potential. For a general cosmological 
fluid, such as the dark matter or the baryons, rj — (3/2a)(<I> + ^ -|- w), where ip is the local gravitational potential (up 
to an appropriate rescaling, whose precise form will be given in Section 2), which must be consistently determined 
via the Local Poisson equation, and w is the specific enthalpy of the fluid (up to the same rescaling), which vanishes 
in the DM case. Therefore, our potential rj in the RHS of eq. (^ is far from being 'external', as it is, for instance, 
dynamically related to the velocity field which appears on the LHS of the same equation. 

The main idea of this paper is that, if?; is given an approximate expression, e.g. by using the results of perturbation 
theory, it can be legitimately treated as a truly external random potential. The choice of variables in eq. (0) is such 
that, if this approximation technique is applied to the DM evolution, to first order it reduces to the well-known 
adhesion model (Gurbatov, Saichev & Shandarin, 1985, 1989; Kofman & Shandarin 1988), which was introduced in 
cosmology to extend the validity of the Zel'dovich approximation (Zel'dovich 1970) beyond the epoch of first caustic 
formation. In the DM case, therefore, only a second-order calculation would produce a non-vanishing external force, 
whose presence would then affect rather small scales. Quite different is the case of the coUisional baryon component, 
where, already at first order in perturbation theory (both Eulerian and Lagrangian) a non-zero contribution to rj 
generally comes out: it is originated from the unbalanced composition of Hubble drag, local gravity and gas pressure, 
thus representing a genuine baryonic feature. 

The presence of the kinematical viscosity term requires some explanation. Its role in the present context is 
twofold. First, it prevents the formation of multi-streams, which are well-known to affect the DM dynamics, but may 
also appear as a spurious effect in the coUisional case, when pressure gradients are given an approximate form in 
terms of linear theory. Second, it allows to transform the problem into a linear one, through the so-called Hopf-Cole 
substitution (e.g. Burgers 1974) u — — 2i/VxlnW, where the 'expotential' U obeys the random heat, linear diffusion, 
equation 

^^ = .V^W(x,.) + !^W(x,.), (2) 
or 2u 

whose solution is expressible in terms of path-integrals (e.g. Feynman & Hibbs 1965). As we will see in Sections 5 
and 6, a suitable approximation technique, valid in the small u case, allows to give the velocity field a simple and 
finite form, more convenient for practical applications. 

The forced Burgers equation has been used to describe a variety of different physical problems, ranging from 
interfacial growth in condensed matter physics, where it is known as the KPZ model (Kardar, Parisi & Zhang 1986) , 
to fully developed turbulence (e.g. Bouchaud, Mezard & Parisi 1995). Later in the paper we will come back to this 
interesting connection and discuss both the analogies and the peculiarities of the cosmological application of this 
equation. 

The stochastic adhesion approximation provides an analytical description of the generation, and subsequent 
merging of shocks, which give rise to the thin network of filaments and sheets in the IGM spatial distribution. Their 
existence is clearly observed in the spectra of high-redshift QSOs, e.g. through the presence of common absorption 
lines in the Lya forest of multiple QSOs, with lines of sight separated by several comoving Mpc at redshifts z ~ 2 — 4 
(e.g. Ranch 1998 and references therein). An important property of our model is that it allows to draw the skeleton 
of the IGM distribution through a straightforward extension of the geometrical technique applied in the free adhesion 
model (e.g. Sahni & Coles 1995 and references therein). The present paper will be mostly devoted to provide the 
physical and mathematical bases for our stochastic adhesion model. Simulations of the IGM large-scale structure will 
be obtained only from the simplified inviscid {i/ = 0) model. In a subsequent paper we will implement our algorithm to 
produce numerical simulations of the IGM distribution and to study the statistical properties of the IGM density and 
velocity fields. Let us stress that approximation schemes like the present one can be particularly useful, as they allow 
to better account for the cosmic variance of large-scale modes, which is poorly probed (especially at low redshifts) 
by the existing hydro-simulations, which are forced to adopt small computational boxes to increase the small-scale 
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resolution [e.g. the discussion in (Viel et al. 2001)]. Even more interesting is the possibility to combine our scheme 
with a hydro-code, using the former to provide the large-scale skeleton of the IGM and the latter to achieve the 
required resolution on small scales. 

The approach most closely related to ours is that recently proposed by Jones (1996, 1999), which aims at 
modelling the nonlinear clustering of the baryonic material. In Jones' model, however, a different set of variables is 
adopted, which docs not allow a direct comparison neither with the present scheme, nor with the Zol'dovich and 
the adhesion approximations, in the coUisionless limit. The most important difference is that, in Jones' model, the 
external random term is identified with the local gravitational potential (p, which is treated as an 'external' one, as it 
is essentially generated by the dominant DM component; moreover, no explicit account for the gas pressure is given. 
In our model instead, the external potential rj is obtained by linearly approximating the composition of Hubble drag, 
local gravity and thermal pressure which act on the IGM fluid elements; in our case r] is an external potential because 
it is determined by a convolution of the initial gravitational potential with the IGM linear filter. 

It should be clear from this introduction that the forced Burgers equation might have wider applications in the 
cosmological structure formation problem. Its relevance (in terms of the closely related random heat equation) in the 
cosmological framework has been first advocated by Zol'dovich and collaborators in the mid eighties (Zel'dovich ct al. 
1985, 1987), as a means to describe the possible origin of intermittency in the matter distribution. The intermittency 
phenomenon consists in the appearance of rare high peaks, where most of the matter is concentrated, separated by 
vast regions of reduced intensity. From the statistical point of view, intermittency in a stochastic process is signalled 
by an anomalous scaling of e.g. structure functions (moments of velocity increments): moments of order p > 2, made 
dimensionless by suitable powers of the second-order moment, grow without bound on small scales (e.g. Gartner J., 
Molchanov 1990; Frisch 1995). This may be viewed as increased non-Gaussianity on small scales, which, in Fourier 
space, appears as a slow decrease in the amplitude of Fourier modes with increasing wavenumber and by a definite 
phase relation between them (Zel'dovich et al. 1985, 1987). Actually, the occurrence of this form of intermittency, 
which seems too extreme and far from our present understanding of the large-scale structure of the Universe, is usually 
obtained under special properties of the noise and only appears at asymptotically late times. Much more interesting 
for the structure formation problem is a second phenomenon, called intermediate intermittency, also described by the 
forced Burgers equation, which consists in the formation of a cellular, or network structure, with "thin channels of 
raised intensity (the rich phase), separating isolated islands of the poor phase" (Zel'dovich et al. 1985). This second 
phenomenon is expected to arise as an intermediate asymptotic situation. 

Can one take advantage of the description of the structure formation process in terms of intermediate inter- 
mittency to predict the specific non-Gaussian statistics which characterizes the nonlinear density field? This is, we 
believe, a challenging issue, which would deserve further analysis. 

The prototype distribution that describes intermittency is the Lognormal one, which naturally arises in multi- 
plicative processes, through the action of the Central Limit Theorem (e.g. Shimizu & Crow 1988). In the DM case. 
Coles & Jones (1991) proposed that a local Lognormal mapping of the linear density field can describe the nonlinear 
evolution of structures in the Universe. Detailed comparison with N-body simulations showed that this model fits very 
well the bulk of the probability density function (PDF) of the mass density field in N-body simulations, for moderate 
values of the rms overdensity a (e.g. Bernardeau & Kofman 1995). Where the Lognormal fails is in reproducing the 
correct PDF for the strong clustering regime (a ^ 1), as well as the high- and low-density tails of the PDF, even 
on mildly nonlinear scales. Moreover, the predicted spatial pattern is too clumpy and poorly populated of extended 
structures to reproduce that of simulations (Coles, Melott & Shandarin 1993). A 'skewed' Lognormal PDF is proposed 
by Colombi (1994), to follow the transition from the weakly to the highly nonlinear regime. 

Quite noticeably, the Lognormal model has also been used in connection with the IGM dynamics. Indeed, the 
semi-analytical model proposed by Bi and collaborators (Bi 1993; Bi et al. 1995; Bi & Davidsen 1997; see also Feng 
& Fang 2000; Roy Choudhury, Padmanabhan & Srianand 2000; Roy Choudhury, Srianand & Padmanabhan 2000; 
Viel et al. 2001), to simulate the low-column density Lya forest, is based on a local Lognormal model, similar to the 
one of Coles & Jones (1991). Comparison with the results of more refined techniques has shown that it provides a 
good fit of the column density distribution in a wide range of values, but it tends to underestimate the abundance 
of lower column density systems (Hui, Gnedin & Zhang 1997). Moreover, the Lognormal model for the IGM tends 
to produce an excess of saturated absorption lines in the simulated transmitted flux, compared with real QSO spectra. 

The plan of the paper is as follows. In Section 2 we review the equations which govern the Newtonian dynamics 
of a two-component fluid of dark matter (DM) and baryons, in the expanding Universe; these are solved in Section 
3 at the Eulerian linear level and under fairly general assumptions on the baryon equation of state. This allows to 
obtain the IGM filter connecting the linear baryon density fiuctuations to the DM ones. A simplified version of our 
model is presented in Section 4: it enables one to follow the combined dark matter and baryon dynamics on weaJily 
nonlinear scales, within the laminar regime. In Section 4 we also give the first-order Lagrangian solution for the 
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baryon dynamics, which is then used to perform numerical simulations of the IGM distribution. In Section 5 we 
discuss the modifications introduced in the baryon dynamics by our improved final model, where we add a kinematic 
viscosity term, to avoid the occurrence of shell-crossing singularities (arising also in the baryon fluid, because of the 
approximate treatment of pressure gradients). This leads to our stochastic adhesion model for the dynamics of the 
intergalactic medium. According to this model, the approximate IGM dynamics is governed by the forced Burgers 
equation for the baryon peculiar velocity field, whoso solution can be expressed as a path-integral. In the physically 
relevant limit of small viscosity, a solution can be found through the standard saddle-point technique. In Section 
6 we discuss how to implement our solution in terms of the first-order Lagrangian particle trajectories previously 
obtained. A geometrical algorithm is also outlined, which allows to draw the skeleton of the IGM distribution, given a 
realization of the gravitational potential and the linear IGM filter. A preliminary analysis of the statistical properties 
of the density field obtained through the stochastic adhesion model is given in Section 7. The concluding Section 8 
contains a brief discussion on the possible applications of our model. 



2 DYNAMICS OF DARK MATTER AND BARYONS IN THE EXPANDING UNIVERSE 

The Newtonian dynamics of a self-gravitating two-component fiuid, made of collisionlcss dark matter and coUisional 
baryonic gas is governed by the continuity, Euler and Poisson equations. The continuity equation for the dark matter 
component (indicated by a subscript DM) reads 

+ 3HpoM + V • (pdmVdm) = , (3) 

where p is the mass density, v = adx/dt the peculiar velocity and H the Hubble parameter at time t; the Euler 
equation reads 

^ ^ ' + (VDM • V) VDM = - V0 , (4) 

where is the peculiar gravitational potential. For the baryon fluid (indicated by a subscript 6), we have 

^ + 3Hpb + V • (pbVb) = , (5) 

and 

3(avb) , , Vpb 

^_ + (vb.V)vb = -V<A- — , (6) 

where p is the pressure. The peculiar gravitational potential obeys the Poisson equation V^(^ — A-kGo^&p, where the 
mass-density fluctuation 5p takes contribution both from dark matter and baryons. Let then /dm and /b = 1 — /dm 
be the mean mass fraction of these two types of matter. We have 

= ^i?ofi()m (/dm5dm + fb5h) , (7) 

where Hq is the Hubble constant, Qom is the closure density of non-relativistic matter (both dark matter and baryons) 
today, (5dm and 5h are the fractional DM and baryon overdensities. In the analytical calculations which follow we will 
neglect, for simplicity, the self-gravity of the baryons, i.e. we will put /dm = 1- 

To close the system we need the IGM equation of state, which can be taken of the polytropic form 

PbfcsT Pb/csro(l-|-5b)^ 

Pb = = — , (o) 

where fcs is Boltzman's constant, 7 the adiabatic index, /i the mean molecular weight (for fully ionized gas with 
primordial abundances it is about 0.6), nip the proton mass and pj, the baryon mean density. In writing the pressure 
term we have assumed the power-law temperature-density relation T = To(l -|- 5by~^ where To = To{z) is the IGM 
temperature at mean density, at redshift z (e.g. Hui & Gnedin 1997; Schaye et al. 1999; McDonald et al. 2000). This 
is adequate for low to moderate baryon overdensity (5h < 10), where the temperature is locally determined by the 
interplay between photoheating by the UV background and adiabatic cooling due to the Universe expansion. The 
underlying assumption for this 'equation of state' is a tight local relation between the temperature and the baryon 
density, which is only true for unshocked gas (e.g. Efstathiou et al. 2000). 

In what follows it will prove convenient to change time variable from the cosmic time t to the scale-factor a 
(e.g. Shandarin & Zel'dovich 1989; Matarrese et al. 1992; Sahni & Coles 1995). This defines new peculiar velocity 
fields u = dtx./da = w/{a^H). Let us also introduce the dimensionless comoving densities p/p = 1 -\- 5 and a scaled 
gravitational potential <p = 2((>/{3a^H^). For the redshift range of interest here (2 w 2 — 4), one can write a^H^ w 
HQQom- A more exact treatment of the DM component would require the use of the growing mode of linear density 
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perturbations, as time variable (e.g. Gurbatov et al. 1989; Catalan et al. 1995); once again, for the range of 

redshifts of interest here one can safely assume D+{t) oc a{t). In the Einstein-de Sitter case the present treatment 
becomes exact. 

We have the following set of equations for the DM component 

DSoM 



Dudm 3 , , v7 ^ 
-^ = -2^(uDM+V^) 



Da 



= -(1 + 5dm)V • udm 



(9) 



where D/Da = d/da + u ■ V denotes the convective, or Lagrangian, derivative w.r.t. our new time variable a{t). For 
the baryons, we have 



Da 

DSb 
Da 



3_ 

'2a 



f 2'ykBTo 



(l + 5b)^~'V5b 



-(l + (5b)V-Ub 



The Poisson equation also gets simplified: 
V (p = . 



(10) 

(11) 

(12) 



3 LINEAR THEORY 

3.1 Dark matter in the linear regime 

Let us start by writing the above set of equations for the DM component in the (Eulerian) linear approximation. The 
continuity and Euler equations, respectively, simplify to 

3 

(5dm = — V ■ Udm ; udm = — — (udm + Vi/?) , (13) 

where, from now on, dots will denote partial differentiation w.r.t. the scale-factor a. The solutions are well known, 
and we will simply report them here. Keeping only the growing mode terms we have 

(5dm(x, a) = a5o(x) ; ip{x, a) = ipo{x) ; udm(x, a) = -V(y3o(x) , (14) 
where So = V^yo. 

3.2 Baryons in the linear regime 

In order to make a similar analysis for the baryons we need to specify the time (or redshift) dependence of the baryon 
mean temperature To, which generally depends upon the thermal history of the IGM, as well as on the spectral shape 
of the UV background.^ We will here consider a simple, but fairly general law of the type To{z) oc (1 + z)" which 
will allow to get exact solutions of the linearized baryon equations. 

At high redshifts, before decoupling, the mean baryon temperature drops like To oc (1 + z); when Compton 
scattering becomes inefficient adiabatic cooling of the baryons implies To oc (1 + z)'^, which makes it practically 
vanish before reionization [this will justify our initial conditions for the evolution of the baryon overdensity in eqs. 
(^^]. As the Universe reionizes, the IGM temperature rises and a different redshift dependence takes place. Various 
types of dependence have been considered in the literature. The linear law To oc (l + z) is often assumed for simplicity. 
According to Miralda-Escude & Rees (1994) and Hui & Gnedin (1997), long after hydrogen reionization has occurred 
the diffuse IGM settles into an asymptotic state where adiabatic cooling is balanced by photoheating, leading to the 
power-law To oc (1 -I- z)", with a = 1/1.76. Much steeper an exponent, a — 1.7, has been recently adopted by Bryan 
& Machacek (2000) and Machacek et al. (2000). More complex is the picture which emerges from hydrodynamical 
simulations of the IGM, where the mean gas temperature appears to retain some memory of when and how it was 
reionized (e.g. Schaye et al. 2000). Observational constraints on To(z) should also be taken into account (e.g. Shaye 
et al. 1999, 2000; Bryan & Machacek 2000; Ricotti et al. 2000; McDonald et al. 2000). 

In view of this variety of assumptions and results it seems reasonable to look for solutions of our equations 
assuming a general exponent a (although, in practical applications we will focus on cases with a < 1). Moreover, 
as we will see later, our model can be straightforwardly extended to any redshift dependence of the mean IGM 
temperature. 

* Similar reasonings would actually also apply to the adiabatic index 7, which we will, however, approximate as being constant 
in what follows. 
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Solutions of the hydrodynamical equations for general values of a have been obtained by Bi, Borner and Chu 
(1992). Many authors (Peebles 1984; Soloveva & Starobinskii 1985; Shapiro, Giroux & Babul 1994; Nusser 2000) have 
given solutions for the case q = 1. In particular, Nusser (2000) has obtained the linear IGM overdensity in the case 
a = 1, extending his calculation to the case of non-negligible baryon fraction (i.e. for < /b < 1), i.e. accounting 
for the baryon self-gravity. We will give here an extensive presentation of this problem for general values of a, both 
because we are going to use the linear baryon overdensity in our nonlinear model for the IGM dynamics, and because 
of the specific form taken by the solutions for our set of initial conditions, which had not been previously obtained 
in the literature. The particular case a = 1 needs to be studied separately; this will be done in the next subsection. 

Let us start by writing the linearized continuity and Euler equations for the baryonic component, 

3 / 2'ykBTo 

— Ub + Vv3 + -—2-- 

2a \ ShQilomfJ-mp 



5b = -V-Ub ; Ub--7^(ub + V.^+ .!L ^'^b) ■ (15) 



Combining these equations together, using the linear solutions for the DM and our temperature-redshift relation 
we obtain, in Fourier space, 

3 • 3 fc^ 3 

(5b(k,a) + — (5b(k,a) -I- ^-;^^S^(k,a) = 2^'5DM(k) . (16) 

The redshift-independent wavenumber kj is related to the comoving Jeans wavenumber kj through 

= a"''Pj = 3//g»o^Mmp(l + .) ^^^^ 
-yksTo 

Only for Q = 1 the two wavenumbers coincide and the Jeans length becomes a constant. 

Equation (|l^) will be solved with the initial, or, more precisely, matching conditions at a = a^, 

c n \ i: n \ k r^ \ k l^ \ SoMi'kjai) 5DA/(k, a) 

flb(k, Oi) = dDM(k, ai) ; db(k, a^) = dDM(k, a,) = = , (18) 

at a 

which are appropriate if the IGM undergoes sudden reionization at Zi (Nusser 2000). 



3.2.1 Case To{z) (X {I + z) 

In the case a = 1, we try a solution of the type — Aa" + Ba, where n and B can be found by substitution back 
into equation (p^. This gives 

X n \ -1/4 r;i Q/4 I A -Q/41 , (5DM(k,a) 

5b(k,a) = a [A^a^^ + w/ j ^ _____ ^ (ig) 



where Q — -^1 — 24fc2/fcj and Ai and A2 are integration constants. From the latter expression we see that on 
large scales the homogeneous part only contains decaying modes; on smaller scales, instead, the homogeneous part 
is characterized by an oscillatory behavior with decaying amplitude. Asymptotically in time one recovers the well 
known solution (e.g. Peebles 1993) (5b(k,a) ~ 5DM(k,a)/(l + fc^/fcj). 

In spite of the absence of growing modes, the homogeneous part of our solution does not necessarily become 
negligible at the time scales relevant to our problem. It is therefore important, to exactly evaluate the homogeneous 
part, by relating it to our initial conditions above. For the two constants of integration we find Ai = _Bi5DM(k, a^) 
and A2 = i325DM(k, Oi), where 

5 + Q f 1 \ f k\'' ^ _ 5-Q f 1 \fk'^'' 

2Q \l + kykj^ 

The baryon peculiar velocity immediately follows from the linear continuity equation. We obtain 



Ub 



\kj) [aj \ 2Q \aj 2Q [aj J 



"P'^ (21) 

i+kyky ^ ' 



3.2.2 General case: Tq{z) oc (1 + z)" 

We will now look for the solution of eq. (|l^), for arbitrary values of a, excluding the case a = 1. In order to find the 
full solution to the inhomogeneous equation (|l6| ), we first notice that its homogeneous counterpart is solved in terms 
of Bessel functions, namely, (5hom(k,a) = y~°'[AiJ'ei{y) + A2j'-&{y)], where a — 1/[2(1 — a)] and we introduced the 
independent variable y = \f2Aa.{klkj')a^^'^°' — ^/2Aak/kj . The full solution of eq. (|lq ) for any a, obtained by the 
standard Green's method, reads 
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5b(k, a) = [AiJ&iy) + A2j^&{y)] + [24a2y~""S55-i,a(j/)] <5DM(k, a) , (22) 

where ^^.^(j/) is the Lommel function. 

By imposing the initial conditions ( p^ ) on our solution, we obtain Ai = BiS-dmO^jCh) and A2 = i325DM(k, a^), 
where 



Bo = 



2sin(7r3) 

l + a 



2sin(7ra) 



vr 



(a+1) 



4a 



1 - 12a(25 - l)2/,i-^"S5a-2,a+i K7-a 



Ja+i + — ( 1 - 125(25 - ■''''5'5a-2,a+i ) J'a 



(23) 



and the Bessel and Lommel functions are evaluated at y = yi — \/^a{k/kj)ay^" . 

Once again, by replacing this expression in the continuity equation, we obtain the linear peculiar velocity of the 
baryons, 

4a "I 

(Bi Ja+i - B2j^i&+i)) + 125(25 - 1)553-2,3+1 y'-^°UDM(k) . (24) 
4a ^ ' 



Ub(k,a) 



3.3 IGM linear filter 

The results of the previous section indicate that one can always express the linear baryon density field as a linear 
convolution of the linear DM one. In Fourier space, one can therefore define an IGM linear filter, Wb(fc,a), such that 

5b(k, a) = WUk, a)5DM(k, a) . (25) 

For the case a = 1 we obtain 



W^.(fc,a) = 0+(^Vf-^"'' 



(si ) 



2Q \aj 2Q \a.. 



i + p/fc; 



2 > (26) 



J 



which corresponds to eq. (A7) of (Gnedin & Hui 1998). Let us look at the main features of the IGM filter, at a 
finite time after reionization. It is immediate to check that on large scales Wh tends to unity: on scales much larger 
than the baryon Jeans length, the baryon distribution always traces the DM, due to the negligible role of pressure 
gradients. On scales much smaller than the Jeans length, instead, the baryon overdensity undergoes rapid oscillations. 
As noticed by previous authors (e.g. Gnedin & Hui 1998; Nusser 2000), and contrary to naive expectations, Wh{k, a) 
has no power-law asymptote on small scales. 
In the general case a 7^ 1, instead we have 

Wb{k,a) = 245^2/-''' S'5a-i,a(y) 

^Zllva) +:r_a(t/)Ja+i(yO) (l ^ 24a2y-^^55a-i.a(yO)] 

^) [(^-a(y) Ja(yO - My)JMy^)) (l - 12«(25 - l)2/,'-'"55a-2.a+i(y0)] • (27) 

yi J sm(7ra) '- ^ ' ■' 

It is worth mentioning that for half-integer values of 5 the Bessel and the Lommel functions in the above equation 
can be written in terms of trigonometric functions. 

Analogous relations hold for the baryon peculiar velocity in terms of the DM one [see eqs. ( pl| ) and ( p^ ) above]. 
Also in this case, of course, Wb ^ 1 in the large-scale limit, as it would be easy to see. The oscillating behavior of Sh 
on small scales is in general more complicated and depends upon the value of a. 

Plots of the IGM linear filter for some values of a are given in Figure 1, at various redshifts. Hydrogen reionization 
is assumed to have occurred at a redshift Zi = 7. Note that, on scales smaller than the baryon Jeans length, acoustic 
oscillations, with decaying amplitude, persist even long after reionization. Note also that different equations of state 
lead to a different small-scale behaviour even if the IGM density is plotted vs. the scaled wavenumber k = k/kj. This 
is because eq. (|l^, which governs the linear evolution of baryon density fluctuations, retains a dependence on the 
slope a in the last term on the LHS, which describes the effect of pressure gradients. 



4 MODELLING THE WEAKLY NONLINEAR IGM DYNAMICS 

With our choice of time variable, the Euler equation, both for the DM and the baryons, takes the simple form 
^ = -V, , (28) 
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Figure 1. Plots of the ratio of the baryon to dark matter density, Wh = <5b/'5DM, versus k = k/kj at different redshifts, for 
various values of a in the mean IGM temperature-redshift relation, To{z) oc (1 -|- z)". Hydrogen reionization is taken at Zi = 7. 



where the potential rj takes a different expression for the two fluids. As the force acting on both fluids is conservative, 
in the absence of initial vorticity the flow remains irrotational (actually, this is only true prior to shell-crossing, for 
the DM component), and we can express the peculiar velocity in terms of a velocity potential, u = V^l?. The potential 
■q in the RIfS of the Euler equation reads 



?7dm = — ($DM + 



(29) 



for the DM and IGM components respectively. 

The dynamical model we introduce here is largely inspired by the Zel'dovich approximation (Zel'dovich 1970), 
and is based on replacing the potential rj, which should be consistently calculated using the full set of equations for the 
two fluids, by a mock external potential, obtained by evaluating its expression within linear (Eulerian) perturbation 
theory. Indeed, this model might be seen as only the first step of a more general, iterative approximation scheme. 

In order to evaluate to linear order the quantity r] one can either compute to first order any single contribution or 
compute its Laplacian by taking the divergence of the LHS of the linearized Euler equation. This yields S/^r; « —V • u, 
and, using the linearized continuity equation, V^rj ~ S . 

We thus conclude that, to first order in perturbation theory the potential V^rj is given by the second 'time'- 
derivative of the linear overdensity. This fact immediately implies that, to first order, ?7dm = (neglecting the 
contribution from decaying modes), while 

7jk{k,a)= [aWi,{k,a) + 2Wb{k,a)] ipo{k) = £{k,a)ipo{k) ■ (30) 
In particular, for a = 1 only the homogeneous part of the 5b solution contributes, and we obtain 



£{k, a) — 



-5/4 



(Q-1)(Q^5) 
16a 



(fli ) 



Q/4 _ 
Bi 



(Q + l)(Q + 5 ) fa 
16a 



-Q/4 . 

B2 



(31) 



with Q = \/i — 24fc2/fcj and the integration constants Bi and B2 as in eq. ( |20[ ) . For a 7^ 1 we get 



£{k,a) = V245^ 
kj 



16q2 



Bi [(45 - 2)J-a+i(2/) + yJ&Mv)] + B2 [(2 - 4q) J-_(a+i)(y) + yj^fa+a) (y)] 

(32) 



245^ j/i-9"(i2a _ 6) [2{& - l)y5'5a-3,a+2 + (1 - 2fi)55a-2,a+i] 
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Figure 2. Plots of the ratio £ = rj/ip versus k = k/kj at different redshifts, for various values of a. The reionization redshift 
is Zi = 7. 



where Bi and B2 are given by eq. (|2^). Plots of the function £{k,a) are given in Figure 2, at different redshifts after 
reionization, for various values of a. 

Therefore, our model is described by the two Euler equations: 

DUOM 



Da 







(33) 



for the DM component, and 

^ = -Vr;. , (34) 
Da 

for the baryons. 

The solution of the above DM equation of motion is well known: it corresponds to the Zel'dovich approximation, 
according to which mass elements move along straight lines, with constant 'velocity' impressed by local linear fluc- 
tuations of the gravitational force at their initial Lagrangian location q: UDM(x(q, a), a) — — Vq<^o(q)- We therefore 
have 



XDM(q, a) = q- aVqiy3o(q) 



(35) 



The trajectories of the IGM fluid elements are instead more complex, as our equation implies a non-zero force 
acting on them, 



pa pT 

Xb(q, a) = q - aVq(/3o(q) - / dr / dr' Vx?7(x(q, r'), r') 

Jo Jo 



(36) 



where we have formally extended the time integration from to a, as, before reionization, when the Jeans length is 
negligible, ?7b(x, a < a^) — 0. 

It will also prove convenient in the following to transform the baryon Euler equation into the following 
Bernoulli, or Hamilton- Jacobi, equation for the velocity potential: 



(37) 



The extreme simplicity of our scheme is shown by the fact that the only information needed to evolve the baryon 
distribution in the weakly nonlinear regime is the IGM filter Wh- It is immediate to realize that, because of our 
derivation of 77b, its expression, V^rj « 5, in terms of the linear baryon overdensity is fully general, i.e. it would apply 
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to general reionization histories [i.e. to general To{z) relations], general baryon equations of state and to the case in 
which the baryon self-gravity is properly taken into account (i.e. the general case < /b < 1). Note that the IGM 
linear filter is simply the ratio of the baryon to DM transfer function at the given time VKb(fc, a) = Ti,{k, a) /ToM{k, a) 
(provided the reionization process has been taken into account in evaluating the baryon transfer function Tb). 

One might think that there is some degree of arbitrariness in the particular choice we made of which terms to 
linearize [those on the RHS of eq. (^s])] and which ones to treat exactly (those on its LHS). The strongest motivation 
for such a choice is its analogy with the procedure leading to the Zel'dovich approximation for the DM component. In 
this sense, the choice is unique; as we will see in the next subsection, in fact, the external force Vi^b obtained above 
and its effect on the IGM trajectories are closely connected to the results of first-order Lagrangian perturbation 
theory: no other choice would have provided such a connection. It should also be stressed that other successful 
approximation schemes for the evolution of the DM component, such as the 'frozen flow' (Matarrese et al. 1992) and 
'frozen potential' (Brainerd, Scherrer & Villumsen 1993; Bagla & Padmanabhan 1994) ones, are indeed based on the 
same choice. The physical reason which makes a linear theory evaluation of rj reasonably accurate is that this quantity, 
for both fluids, contains terms that receive their dominant contribution from small wavenumbers. This was indeed 
the original motivation which led Zel'dovich to obtain his celebrated algorithm, although in modern language the 
Zel'dovich approximation is more commonly explained within the first-order Lagrangian approximation (e.g. Sahni 
& Coles 1995 and references therein). The link between these two alternative derivations is discussed in the next 
subsection. 



4.1 Lagrangian approximation to the baryon trajectories 

The trajectory of any mass element can be written in the general form 

x(q,a) =q+S(q,a) , (38) 

where S is the 'displacement vector'. As far as the evolution of the fiuid is far from the strongly nonlinear regime, 
the displacement vector can be considered to be small, i.e. the Eulerian and Lagrangian positions of each particle 
are never too far apart. This consideration is at the basis of Lagrangian approximation methods. The basics of the 
flrst-order Lagrangian scheme applied to our two-component fluid are reported in Appendix A. Applying these ideas 
to our external force, we obtain 

Vxr?b(x(q,a),a) = VqJ7b(q,a)+0(s2) . (39) 

To flrst order in the displacement vector we can therefore replace the Eulerian force Vx?7b(x, a) by its Lagrangian 
counterpart Vq77b(q, a) in the baryon trajectories, which leads to the much simpler form 

Xb(q, a) = q- aVq?/>b(q,a) , (40) 
where the 'baryon potential' tph is defined by 

V^^b(q,a) = Mai^ (41) 

and is related to the peculiar gravitational potential by the Fourier-space expression ■i/;b(k,a) — Wh{k,a)ifoO<-)- 

Moreover, as shown in Appendix A, the trajectories described by eq. (^o|) represent the result of flrst-order 
Lagrangian perturbation theory applied to our baryon gas. Similar results have been recently obtained by Adler & 
Buchert (1999) for the case of a single self-gravitating coUisional fluid (i.e. for /b = 1). There is an important difference 
between these trajectories and the DM ones. According to the Zel'dovich approximation, DM particles move along 
straight lines with constant velocity, whereas the baryons are generally accelerated along curved paths; this is due to 
the non-zero force, resulting from the composition of three terms: the local Hubble drag, the local gravitational force 
caused by the dominant DM component and the gradient of the gas pressure. 

The baryon peculiar velocity which follows from this first-order Lagrangian approximation is 

Ub(x(q,a),a) = -Vq (aV>b(q,a) +7/)b(q, a)) , (42) 

which, unlike the DM one, deviates from the initial velocity Ubo(q) = Ub(x(q, 0),0) — — Vqiy5o(q). 

Let us mention an important property of the approximate expression for the velocity field that we obtained 
using flrst-order Lagrangian perturbation theory: unlike what happens in the DM case, the vector Ub of eq. (^2|) is 
irrotational in Eulerian space only if we restrict its validity to flrst-order in the displacement vector; a non-zero velocity 
curl component, Vx x Ub 7^ 0, in fact arises beyond this limit. This feature is not expected to imply serious problems 
as long as the system has not evolved too deeply into the nonlinear regime, i.e. until the baryon overdensity has not 
reached values ^ 1. This is an unphysical feature deriving from the extrapolation of the first-order results to a regime 



© 0000 RAS, MNRAS 000, 000-000 



The growth of structure in the intergalactic medium 11 



where higher-order terms should be taken into account. A similar feature appears in the Lagrangian perturbation 
approach to the DM dynamics, when both the growing and decaying solutions are included (e.g. Buchert 1992). No 
problems of this type, however, occur when the Eulerian scheme above is adopted, i.e. when eq. (^) and its solution 
( p6[ ) are assumed. For this reason we will consider our 'Eulerian model' of eqs. (|3^) and ( p6| ) as the correct one and 
the Lagrangian scheme leading to eq. ( po[ ) as being essentially a shortcut to get approximated baryon trajectories. 

The slight discrepancy between the Eulerian and Lagrangian schemes used to derive the external force Vrj acting 
on the baryons is a peculiarity of the coUisional case. In the DM case, these two techniques - linearizing the RHS of 
eq. ( p^ in Eulerian space, and expanding to first order in Lagrangian perturbation theory - lead to identical results. 

Approximation schemes for the low-density IGM dynamics which are closely related to our Lagrangian treatment 
have been studied by Reisenegger & Miralda-Escude (1995), Gnedin & Hui (1996), and Hui, Gnedin & Zhang (1997). 
There are, however, important differences that we would like to point out: i) in our model the gas trajectories, and 
thus the resulting IGM spatial clustering depend on the specific ionization history (different values of a in the simplest 
case); ii) our model is by definition able to exactly reproduce the behavior of the baryon component in the linear 
regime. In practice, while previous models adopt an IGM filter which is best-fitted to simulations, ours directly derives 
from baryon dynamics. 



4-1.1 Numerical simulations of the IGM distribution in the laminar regime 

In order to display the effect of the IGM filter on baryon trajectories, we have produced a set of numerical simulations 
based on the Lagrangian scheme. 256^ particles were laid down on a uniform cubic grid and moved according to 
eq. (^o[). Realizations of the peculiar gravitational potential were obtained in Fourier-space according to standard 
algorithms. The model shown in Figure 3 is a spatially flat, vacuum dominated cold dark matter (CDM) model, 
with h = 0.65 (where h is the Hubble constant. Ho, in units of 100 km s~^ Mpc~^), a cosmological constant with 
present-day density parameter floA ~ 0.7, CDM with rJoCDM = 0.26 and baryons with the remaining floh = 0.04. 
A Zel'dovich primordial power-spectrum of adiabatic perturbations is assumed, with the standard normalization of 
ACDM modelsf], as = 0.99 (Viana & Liddle 1999), where ag is the rms mass fluctuation in a sharp-edged sphere of 
radius 8 Mpc; we adopt the Bardeen et al. (1986) CDM transfer function, corrected to account for the baryon 
contribution, as in (Sugiyama 1995). 

The box-size is 10.24 Mpc. Three different IGM thermal histories are considered; in all cases we assume that 
reionization occurred at Zi = 7, but we assume that the mean IGM temperature evolves according to different power- 
laws afterwards, To{z) oc (1 -|- z)", with a = 1.0, 0.8, 0.6. Figure 3 shows a thin (0.02 Mpc) slice of our simulation 
box ai z — 3. The corresponding DM distribution is also shown for comparison. Note that different values of a lead 
to a different small-scale distribution of the IGM, because of the different time dependence of the Jeans scale and 
ii) the different trajectories followed by the baryons after reionization. The comoving Jeans wavenumber kj has been 
set to the same value (7 Mpc~^) in all our models, at redshift z — 3. 



4-1.2 The IGM density field in the laminar regime 

Once particle trajectories are known, it is immediate to obtain the density fleld, using mass conservation, l-|-5b(x, a) = 
||9xb(q, a)/9q||~^. If we adopt the Lagrangian expression of eq. (^o|), we get 

(43) 

The strain tensor d^tpb/dqidqj can be locally diagonalized along principal axes, whose direction will generally 
depend upon time, unlike the DM case. Calling Ai the corresponding eigenvalues, we write 

3 

1 + <5(x(q, a),a)^'[[[l- a\, (q, a)]-\ (44) 
1=1 

which shows that a caustic singularity will form at the finite time asc(q) = 1/A*(q, asc), where A*(q, asc) is the 
largest eigenvalue of the time-dependent strain tensor. The time dependence of the eigenvalues, which on small scales 
becomes oscillatory, implies that the shell-crossing condition can even be met more than once by a given mass element 
along the same principal axis. 

t To correct for the small inaccuracy introduced by the use of a, instead of the DM growth-factor, as time variable in our 
treatment, we apply a correction factor (1 -|- z)D+{z) to the 0"g normalization, where D+{z) is the linear growth-factor of DM 
density fluctuations, normalized to unity at 2 = 0. 



1 + (5b(x(q, a),a) 



det 



^ i9^V'b(q, g) 
dqidqj 
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Figure 3. The spatial distribution of the dark matter and IGM, with different temperature-redshift relations, are shown at 
redshift z = 3. All models have the same Jeans length, kj = TMpc"-"^, at this redshift. 



The extrapolation of our Lagrangian approximation beyond the formation of the first pancakes, leads to an 
artificial diffusion of these structures. This problem becomes more and more severe at low redshift, making any 
simplified description of the IGM - and of the Lya forest - based on the Lagrangian trajectories quite unreliable. 
The stochastic adhesion model presented in the next section aims precisely at overcoming this problem. 

One might use eq. ( ^ ) for the density field to obtain an analytical expression for the probability distribution 
function (PDF) of the IGM density field, by a simple extension of the traditional Doroshkevich (1970) formalism 
(e.g. Kofman et al. 1994). Such a technique has been applied by Reisenegger & Miralda-Escude (1995), who adopted 
a smoothed Zel'dovich approximation for the IGM evolution. Our model, however, allows to get a more refined 
description in which the PDF is sensitive to the IGM thermal history. Similarly, one might also obtain statistical 
information on the Lya forest in quasar absorption spectra, e.g. in terms of the column density distribution function 
of the lines (e.g. Hui, Gnedin & Zhang 1997); this too would display a dependence on the assumed To{z), which might 
be tested against observational data. These, and related applications of our Lagrangian algorithm will be presented 
elsewhere. 

The emergence of shell-crossing singularities in the velocity pattern of our coUisional fiuid should be understood 
as an artifact of the linearized treatment of pressure gradients in the Euler equation. This feature can be also seen 
as follows. If we evaluate the force rj on large scales, to lowest order in k/kj we find, for q = 1, 



1/2 / k_ 
J . 

which comes only from the homogeneous part of eq. dl9L and 



for Q 7^ 1, from the inhomogeneous term of eq. (^^. In both cases the result can be expressed in the form Vr; ~ 
— //cff (ffl)V^uo, where VcS > 0, for 0.5 < a < 1. This shows that, as long as the velocity is in the linear regime, the 



© 0000 RAS, MNRAS 000, 000-000 



The growth of structure in the intergalactic medium 13 



effect of the linearized pressure gradient is similar to that of an effective kinematical viscosityQ with coefficient Vca, 
which smooths the velocity field. As soon as the system enters the mildly non-linear regime, however, u deviates 
from its initial value and this smoothing is no more effective in stopping the formation of caustic singularities in the 
densest regions. 

Two processes will actually prevent the formation of multi-streams in our physical system: the first is due to the 
binding action of gravity, which tends to keep pancakes and filaments thin, an effect which is experienced both by the 
DM and IGM components; the second is the nonlinear action of the gas pressure, which is instead characteristic of 
the baryonic component. Because of the difficulty to deal analytically with the gas pressure beyond any perturbative 
approximation, we will find convenient in the next section to introduce an artificial viscosity term in the baryon Euler 
equation, whose effect is to smear these shell-crossing singularities of the velocity field. 



5 THE STOCHASTIC ADHESION APPROXIMATION 

The model for the formation of structures in the baryon distribution discussed in the last section is based on a 
suitable approximation for the force exerted on fluid elements by gravity and by the surrounding fluid patches 
through pressure gradients. The latter are most important on small scales where their effect is to smooth the gas 
fluctuations relative to the DM ones, thus preventing the occurrence of shell-crossing singularities in the coUisional 
component; nonlinear effects here manifest themselves through the formation of shock- waves. In our scheme, just like 
in all similar schemes [e.g. the modifled Zel'dovich scheme used by Gnedin & Hui (1996) and Hui, Gnedin & Zhang 
(1997)], the gas pressure is only included through a linear approximation. Because of this fact, when the system enters 
the nonlinear regime, i.e. when two particles come very close in space, shell-crossing can affect their evolution leading 
to the subsequent occurrence of multi-streams. In order to prevent this unphysical phenomenon we will adopt the 
same technique which proved so successful in the DM case to extend the Zel'dovich treatment beyond its actual range 
of validity: we add a kinematical viscosity term to our approximate Euler equation. This method is at the basis of 
the adhesion approximation for the formation of large-scale structure in the coUisionless case (Gurbatov et al. 1985, 
1989; Kofman & Shandarin 1988). The adhesion model is based on the three-dimensional generalization of Burgers' 
equation of strong turbulence. It is the simplest equation which describes the formation and subsequent merging 
of shocks. According to the adhesion model, DM particles move according to the Zel'dovich approximation until 
they fall into pancakes, when, owing to the viscous force, they stick together. Next, pancakes drain into fllaments, 
and fllaments into clumps. The thickness of these structures is monitored by the value of the kinematical viscosity 
coefficient v, and becomes infinitely thin as u ^ 0. 

We assume that the equation of motion which governs the IGM dynamics is (from now on we will avoid the 
subscript 'b' on baryon quantities, where unnecessary) 

£^ = u\/'u-\/v, (47) 
Da 

where the kinematical viscosity coefficient v is here assumed to be small, but non-zero. In our coUisional case, moreover, 
there can be an extra reason to add such a term: the gas is effectively experiencing some shear viscosity, although on 
scales much smaller than those under consideration. The physical viscosity term should actually depend upon space 
and time, through the inverse of the local gas density. One might however argue that, in the physically relevant limit, 
where only a tiny viscous force is present, even a constant i/ will produce the correct qualitative trend. It might be 
interesting to mention an alternative interpretation of the viscosity term in self-gravitating systems. According to 
Dominguez (2000), a term which resembles the one for kinematical viscosity is the unavoidable consequence of the 
coarse-graining process inherent in a hydrodynamical description. 

The equation above is known as the forced Burgers equation of nonlinear diffusion (e.g. Barabasi & Stanley 1995; 
Frisch & Bee 2000). Its possible application to the cosmological structure formation problem was suggested long ago 
by Zel'dovich et al. (1985, 1987) and recently studied in greater detail by Jones (1996, 1999). There are a number 
of differences between our approach and the one by Jones, which lead to a different dynamics of the IGM. First, 
we used a different time variable (the scale-factor a instead of the cosmic time t) to define the baryon velocity field 
and its acceleration, thus leading to the form of eq. (^^, where the force on the RHS was replaced by its linearized 
expression. Second, we do include some effect of the baryonic pressure in our external random potential 77, which is 
instead absent in Jones' treatment, where the role of rj is played by the linear gravitational potential generated by 
the dark matter component. Thus, our scheme, unlike the one by Jones, is able to exactly reproduce the evolution 

t A related point has been recently made by Buchert & Dominguez (1998). 
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Figure 4. The evolution of the ratio £ = rj/ipo, for various models of the IGM mean temperature evolution. 



of the baryons at the Unear level. Finally, our random potential has a non-trivial time dependence, unlike the one 
assumed by Jones. 

The forced Burgers equation can be transformed into a Bernoulli-like equation for the velocity potential, namely 
9$ 1 



9a + 2(^*)' = 



V 



(48) 



The problem is fully specified once the initial conditions for the velocity potential and the statistics of the noise 
are given. Our initial velocity potential is "l'o(q) ~ — ¥'o(q)- The statistics of the noise follows directly from that of 
the linear gravitational potential, which we assume to be a Gaussian random field. Then, also is a Gaussian process 
with zero mean and auto-correlation function 



{ri{x,a)i]{x + r,a')) 



1 

2^ 



dk 

fc2 



£{k,a)£{k,a')P{k)jo{kr) , 



(49) 



where the average is over the ensemble, -P(fe) is the linear power-spectrum of DM density fiuctuations, extrapolated 
to the present time, and jo is the spherical Bessel function of order zero. According to this formula, 77 is a stochastic 
process with the following properties: 

It is a colored (i.e. non- white) Gaussian random process both in space and time. 
a) It is stationary (i.e. homogeneous and isotropic) in space but not in time (as its auto-correlation function does 
depend on |x — x'| only, but not on |a — a'|: this is typical of the cosmological case. 

Hi) On small scales, 77 oscillates rapidly in time for all a < 1, with a period which generally depends on the wavenumber. 
iv) The non-separability of the space and time dependence of its auto-correlation function implies that 77 behaves as a 
stochastic process both in space and time. More precisely, the time evolution of the noise in each given spatial point 
cannot be predicted on the basis of local initial conditions only, as it depends on the realization of the Gaussian field 
ipd over the whole Lagrangian space, through the time evolution of the IGM filter. 



Let us finally stress an important peculiarity of our cosmological application of the forced Burgers equation: 
the same stochastic process ^0 which underlies the external random potential rj — £ * ipo also provides the initial 
condition for the velocity potential, $0 = —<Po- 

A plot of the function £{k, a) vs. redshift is given in Figure 4 for various values of a, to display the time-dependence 
of the random force and, in particular, its oscillatory character on small scales. In Figure 5 the dimensionless power- 
spectra A^(fc) of two terms contributing to the evolution of the velocity potential are shown: that of the random 
potential rj, A'^{k,z) = {l/2Tv^)k~^£^{k, z)P{k), and that of the linear velocity potential $0, multiplied by a factor 
l/a = l + 2to give an estimate of the 'deterministic' term in eq. (^), namely Ajj^_^^j^^ {k, z) = {l/2n^){l+z)^k~^ P{k). 
Note that the two contributions become of the same order at fc ~ fcj at the relevant redshifts; on smaller scales the 
random potential dominates the IGM dynamics. 

The above evolution equation for "1>, in cases when the random potential rj is white-noise in time, has been 
extensively studied in condensed matter physics, where it became popular as the Kardar-Parisi-Zhang (KPZ) equation 
(Kardar, Parisi & Zhang 1986). The KPZ equation describes the growth of an interface under random particle 
deposition. Here the viscosity coefficient v plays the role of temperature and the velocity potential $ is interpreted 
as the height above an initially flat surface, which is driven by the random noise and gradually becomes rough. 
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Figure 5. The dimensionless power-spectra of the stochastic and linear deterministic terms in the evolution of the IGM velocity 
potential are shown for o = 1 at z = 3 (left panel), z = 2 (middle panel) and 2 = 1 (right panel). The normalization of the 
linear density power-spectrum, P{k), is chosen arbitrarily. 



5.1 Solving the random heat equation 

Equation ( ji^ ) can be easily related to a linear, partial-differential equation, by means of the nonlinear Hopf-Cole 
transformation $ = —2v\nU (e.g. Burgers 1974), which leads to the linear diffusion equation with a multiplicative 
potential, also called 'random heat' equation 

^ ^ = W(x,a) + ^- — '-U{-x.,a) . (50) 

Starting from the solution of the latter equation for the Hopf-Cole transformed velocity potential U [which has 
been dubbed 'expotential' by Weinberg & Gunn (199Qa)] one can easily find the velocity by transforming everything 
back 

u = -2.^ (51) 

One first obtains an expression for the 'transition kernel' /C(x, a|q, 0), representing the particular solution obtained 
from a Dirac delta function, 5(x — q), at the initial time a = 0. This has a formal solution given in terms of the 
(Euclidean) Feynman-Kac path-integral ^ formula (e.g. Feynman & Hibbs 1965), namely 

/C(x,a|q,0)= I ' \Dx(r)]e-^''^% (52) 

J q 

where the action is given by 

(53) 



S[x] = \ dr/:(x,x,r) ^ dr 
lo Jo 



with £. the Lagrangian of a particle moving in the potential 77. 

In our diffusion equation the transition kernel is immediately understood as the conditional probability of finding 
a particle in x, at time a, given that it was initially in the Lagrangian position q. To better understand the path- 
integral solution it is however convenient to think it in connection with a quantum mechanical problem. It is in fact 
immediate to realize that an inverse 'Wick rotation' of the time variable a ia, together with the formal replacement 
21/ ^ h transforms eq. ( |5o| ) into the Schrodinger equation for a particle, of unit mass, subjected to the potential 77. 
According to the path-integral representation of quantum mechanics, its solution is obtained by 'integrating' over all 
possible paths connecting these two end-points, each one weighed by the action S[x.], calculated along this path. 

Once the kernel is known, the solution of the random heat equation with appropriate initial conditions is 
obtained through the application of the Chapman-Kolmogorov equation (e.g. van Kampen 1992), by integrat- 
ing the product of the initial function Uo with the transition kernel over the whole Lagrangian space, namely 
U{x,a) = J d^q Wo(q) A;^(x, a|q, 0). In our case Wo(q) = exp[— <l?o(q)/2i'] — exp[(y3o(q)/2i^]. Thus, 

W(x,a) = / d^g e~*«<'^'/^'-' / ' ' [I?x(r)] e^®''^'' . (54) 



q 



The solution of eq. (^) in the absence of a potential, i.e. for the 'free' adhesion approximation, is reviewed in Appendix B. 
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In the limit of vanishing viscosity (corresponding to the classical limit, ?i ^ 0, in our quantum mechanical 
analog) the dominant contribution to the path-integral comes from the 'classical' path, i.e. that which satisfies the 
Euler-Lagrange equations of motion, 

dC_d_d£^^^ (55) 
9x da 9x ' 

which in our case leads to Newton's second law, x = —Vrj. Thus, the particle trajectories along the classical path 
read 



pa pT 

Xci(q, a) = q + auo(q) - / dr dr'V 77(x(q, r'), r') 
Jo Jo 



(56) 



with general initial velocity uq. Note that there are still infinitely many classical trajectories joining the two end-points 
(q, 0) and (x, a), corresponding to the freedom to choose the initial velocity uo(q). 

We can now expand around the classical trajectories in the standard manner of quadratic approximations (e.g. 
Feynman & Hibbs 1965), that is x = Xd + subject to the constraint that the fluctuations around the classical 
trajectory vanish at the end points, namely ^(0) = ^(a) = 0. We then have 



with the symbol S standing for functional differentiation. The classical action, ^d, is a function of the end-points 
(q, x) and of the time interval a, which has to be calculated by replacing the solution of the Euler-Lagrange equations 
into its expression. 

The action is an extremum for the classical trajectory, thus (55'/(5x)x^i ~ 0. Inserting the expansion ( [57| ) in the 
path-integral gives 

/C(x,a!q,0)=e-^^'/^'' / [©CW] e ^ ^====01, (58) 







where we used the fact that the Jacobian of the transformation from x to ^ is unity. Note that the path-integral in 
( ^8| ) starts and ends at zero, because the fluctuations vanish at the end-points. Contrary to an ordinary integral, the 
fact that the upper and lower limits are the same does not imply that the path-integral vanishes. This can be clearly 
seen by considering the path-integral of a free particle (see Appendix B), which can be evaluated exactly and does 
not vanish even if one is only interested in the path which starts and ends at the same point (e.g. Feynman & Hibbs 
1965). Since the integral over ^ starts and ends at 0, it can only be a function, F{a), of the end-times. Therefore, we 
can rewrite the above path-integral as /C(x, a|q, 0) — F(a) exp[— 5'ci/2!^]. We need not find the explicit expression for 
the pre-factor F{a), since this function cancels for the velocity. The solution of the random heat equation, for small 
1/, reads 



with 

$ci (x, q, a) = Sci (x, q, a) -I- $0 (q) . (60) 

Replacement into eq. (^l|) yields 

, , Jd'q Vx&i(x,q,a) e-^'o'(--''-°)/^- 
"(''' = Jd^q e~^'c.(x,...)/2. ■ (61) 

In the case of vanishing force, eq. (^l|) reduces to the well-known exact solution of the three-dimensional Burgers 
equation (see Appendix B), which underlies the adhesion approximation for the dark matter component. In our case, 
it represents a useful approximation, which could be used for a numerical evaluation of the baryon velocity field, 
by extending the technique applied by Nusser & Dekel (1990) and Weinberg & Gunn (1990a,b). The main practical 
difficulty in using the solution (^) is, however, represented by the lack of an explicit expression for the force Viy, 
as a function of x and a. A possible shortcut would be to exploit the Lagrangian approximation for the trajectories, 
introduced in Section 4.1. This method will be discussed in Section 6. 

Once the peculiar velocity field of eq. ( |6ll ) is known at each Eulerian point x as a function of time, one can obtain 
the actual baryon trajectories, by numerically integrating the integral equation (e.g. Nusser & Dekel 1990; Weinberg 
& Gunn 1990a,b) 

x(q,a)=q+/ dr u(x(q, r), r) , (62) 
Jo 
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starting from Lagrangian particle positions on a grid. From these trajectories the baryon density can be found, 
according to the standard formula 1 +S{x,a) = ||9x(q, a)/9q||~^. 



5.2 Steepest-descent approximation 



In the limit of small u we can evaluate the integral over initial positions in eq. (|5^) using the steepest-descent, or 
saddle-point, approximation. We have 



W(x,a) = i^(a) e-*-'<"'''='">''2''y d^<5g e 

= F{a) {ATTuf^Y, J. (x,q.,a)e-*=' 

where 

a^$ci(x,q,a) 



js(x,qs,a) = det 



dqidqj 



-1/2 



(63) 



(64) 



The sum in eq. (|6^) extends over the saddle points q^, i.e. the Lagrangian coordinates corresponding to the 
absolute minima of the function $ci(x, q, a), for fixed x and a. These are found by solving the equation 



Vq$ci(x,q,a) 



= Vq(Sci(x,q,a) +$o(q)) 



= 



(65) 



Since, for the classical action, VqSd — — uo(q), we see that the classical orbit ( |56[ ) passes through the saddle point 
qj, if the initial velocity satisfies uo(q) = Vq$o(ci) and $(x, q, a) > <3?(x, q^, a), for any q 7^ q^. 
We then obtain 



x(qs, a) = qj, 4- aVq'l>o(qs) 



dr / dr' Vx??(Xci(qs,T'),r') 



(66) 



for the particle trajectory, which coincides with eq. (Pq), as $0 = — V'o- 

Replacing eq. ( |63| ) in eq. (^) we finally obtain our saddle-point solution of the forced Burgers equation for the 
baryon velocity field 



u(x, a) = ^ Ws(Xiqs,a)VxSci(x, qs,a) , 

S 

where the (normahzed) weights read 

j.(x,q.,a)e-*<=.(-.q..-)/2- 



Wsl^, qs, aj 



js (x,q3,a) 



(67) 



(68) 



(69) 



J^s js (x, qs , a)e-*ci (x,q. ,a)/2. (x, , a) 

(the last step is justified by the fact that, for a given x and a, all the absolute minima have the same "l>ci)- 

Since the Eulerian gradient of the action is the velocity along the classical trajectory, we have from eq. ( p6| 

VxSci(x, qs,a)) = Uci(x(q,a),a) = / dr Vx77(x(qs, r), r) -|- dr dr' Vxr;(x(qs, r'), r') . 

" Jo Jo Jo 

At early times, there is a single contribution to the steepest-descent velocity field, coming from a unique La- 
grangian saddle point q^ that minimizes $ci for given x and a, so that u(x, a) — u(xci(qs, a), a). With time passing 
the mapping q^ ^ x becomes many-to-one. At this stage, according to the steepest-descent solution, the fiow at x 
becomes a weighted average of the velocity of all the mass elements converging to that point at a along their classical 
orbits. This sets the onset of the epoch of pancake formation in our model. This interpretation of our solution will 
become more clear from the geometrical construction outlined in the next Section. 



In concluding this Section, let us come back to our analogy with the KPZ model. A different version of it, which 
is also often used in condensed matter physics, is obtained by making the Hopf-Cole nonlinear transformation on the 
KPZ equation, which leads to the linear diffusion equation for the free energy Z = exp(— $/2i^) (Kardar & Zhang 
1987). Its path- integral representation is then interpreted as follows. The path traced by a particle starting from an 
initial point x(0) = q and arriving at a final point x(a) — x, moving in D spatial dimensions, can be viewed as a 
polymer in D -I- 1 dimensions joining (q, 0) to (x,a), with the time variable interpreted as a spatial one along the 
direction of main extension of the polymer. One specifically uses the term directed polymer to emphasize that the 
path is constrained to go only forward in time without crossing itself. The case studied here of vanishing viscosity is 
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Figure 6. In the top figures the path of a single particle before falling into a shock is shown for the stochastic (IGM) and free 
(DM) adhesion models. The central right panel shows a collection of inertial paths which form when particles collide, merge 
and move together, conserving momentum. The left central figure shows an ensemble of paths (polymers) of lowest energy in 
a random energy landscape. Each polymer has its end points fixed and finds the optimal path in between. The shocks of the 
Burgers equation are the optimal polymers in the language of the interfacial growth. The tree structure is a characteristic of 
many optimization problems. The bottom panels show the velocity field at a fixed large time in the three cases of zero viscosity 
(Zel'dovich), free adhesion and stochastic adhesion. In the Zel'dovich approximation, multivalued regions bounded by caustics 
are formed since the particles can move through each other. The caustics are replaced by shocks of finite mass once a viscosity 
term stops the shell-crossing of the particles. This holds for both free and stochastic adhesion models. Finally, a new feature is 
introduced into the picture: in the presence of a stochstic force, particles trajectories are no longer inertial, but are curved, as 
shown in the bottom left panel. 



frequently referred to, in the language of KPZ and directed polymers, as the 'strong coupling', or 'zero-temperature', 
limit. This coincides with the limit of large Reynolds number, in the language of turbulence, or to the classical limit, 
in our quantum mechanical analog. 



6 DRAWING THE SKELETON OF THE IGM DISTRIBUTION 

There are many possible ways to look for a numerical solution of the forced Burgers equation. Interesting techniques 
might be devised, based on the analogy with the KPZ model. In the v —> limit, in fact, the problem of finding the 
'free energy' Z, corresponding to our potential W, reduces to a simple optimization problem. In the KPZ context, 
one frequently uses a numerical scheme referred to as the transfer matrix method (Kardar & Zhang 1987). First, 
one generates a random energy landscape (e.g. in 2D, on a square lattice, one puts random numbers with a certain 
distribution and correlation, depending on the statistical properties of the noise rj) . The optimal particle starts at a 
lattice point q, chooses the next lattice point which would cost it the least energy, jumps to it and so on, continuing 
its journey to the final point x. In this way an optimal polymer, joining (q, 0) to (x, a), is formed with the minimum 
amount of energy. The free energy Z is then evaluated by weighing the optimal polymer, which ends there, over all 
the possible starting points. A ID illustration of this numerical method is given in Figure 6. 

The application of a similar scheme to the cosmological structure formation problem is however complicated by 
the colored and non-stationary time-correlation properties of our random potential, as well as by the fact that the 
cosmological interesting regime does not usually occur at asymptotically large times, unlike what happens in the case 
of directed polymers. 

Let us next describe an alternative method, based on our steepest-descent solution of the forced Burgers equation. 
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Two main technical difficulties arise, when one tries to implement the saddle-point solution to simulate the baryon 
distribution. The first concerns the velocity along the individual classical trajectories, which contribute to the sum in 
eq. (p^: this involves a step-by-step integration of the force along the trajectory. The second concerns the weight Ws 
which requires explicit knowledge of Sci along these classical trajectories. The problem is equivalent to that of finding 
the velocity potential, $ci(x, q, a) = S'ci(x, q, a) + <l?o(q), which solves the 'classical' Hamilton- Jacobi equation (^) 
with initial conditions <l>ci(q, q, 0) = $o(q). The technique is therefore completely analogous to that used in the free 
adhesion approximation for DM particles, but is made more involved by the presence of the random potential 77. 

A substantial simplification is obtained, if we adopt the same approximation of Section 4.1, which basically 
consists in expanding the particle orbits to first order in the displacement vector. Let us describe the method in 
detail. For each orbit the system contains a cyclic variable corresponding to the initial velocity, namely 

uo = + Vq(V - ^0) , (70) 

where tp is the baryon potential defined in eq. (^l|). Replacing x = uo -f- aVq((/3o — V') in the action (|5^) and then 
using (^^, we obtain, to second order in the displacement vector S, 

fx )^ 1 2 

$,i(x, q, a) = - a(x - q) • Vqi/'(q, a) - 7/;(q, a) - a^(q, «) " 3 _^ '^^^^ (VqV'(q, r)f + 0(8^) , (71) 

where, for consistency, we expanded the potential -q to first order in the displacement vector around the Lagrangian 
coordinate (as 77 itself is a first-order quantity), i.e. 

77(x,a) =77(q,a)+Vq77(q,a)- [auo(q)- Vq(V(q,a)-^o(q))]+0(S') , (72) 

It is immediate to check that the velocity potential of eq. (|7^) gives the required solution of the Hamilton- Jacobi 
equation (p7|), if the potential rj in its RHS is consistently expanded as in eq. (|7^). From eq. ( |7l| ) we can immediately 
obtain an explicit expression for js which appears in our saddle-point solution through the weight of eq. (^). We 
have 

\ -1/2 

Sij 9^?/)(q, a) 9^t/)(q, a) 



js(qs, a) = det 



4- a- 



+ 0(8'). (73) 



dqidqj dqidqj 
The Lagrangian approximation for the classical trajectory 
Xci(qs,a) = qs - aVq7/'(qs,a) (74) 
and the corresponding velocity 

Uci(x, qs, a) = — 3i _ t(Vqi/'(qs,a) (75) 
a 

are the remaining ingredients of the saddle-point solution (|6^). Note that everything is now explicitly given in terms 
of the linear baryon potential and its time derivatives. 

We are now ready to introduce a geometrical construction of the saddle-point solution, obtained in analogy with 
the adhesion approximation, which aims at obtaining the skeleton of the IGM distribution. The graphic technique 
has been largely applied in connection with the solution of the one- dimensional Burgers equation (e.g. Burgers 1974). 
In the cosmological context this method has been described in a number of papers (e.g. Gurbatov et al. 1985, 1989; 
Kofman & Shandarin 1988; Shandarin & Zel'dovich 1989; Vergassola et al. 1994; Sahni & Coles 1995) and applied in 
one, two and three-dimensional numerical simulations of the large-scale structure of the Universe (Kofman, Pogosyan 
& Shandarin 1990; Williams et al. 1991; Kofman et al 1992; Sahni, Sathyaprakash & Shandarin 1994; Sathyaprakash 
et al. 1995), finding good agreement with N-body simulations in models which have sufficient large-scale power (e.g. 
Sathyaprakash et al. 1995). 

Let us start by noting that, for given x and a, the classical velocity potential <l>ci(x, q, a) describes a random 
three-dimensional hypersurface in Lagrangian space. The idea is that absolute minima of this hypersurface correspond 
to Lagrangian points of first-touching with the hyperplane E(i(q) = h, as the parameter h is raised from —00 |^. 

The geometrical construction can proceed along similar lines as for the adhesion model, namely: 
i) At each time a, and for each point x in Eulerian space, find the hyperplane Sh(q) tangential to <l?ci(x,q, a) at 
the Lagrangian point(s) q^^, by increasing the parameter h from —00, so that the two hypersurfaces do not cross 
anywhere. 

a) Once the saddle point(s) q^^ (x, a) are found, the corresponding velocity at x at time a is given by eq. (|67|), as the 

^ In the free adhesion case it is actually more convenient to look for the first-touching points of the paraboloid P/i(q|x, a) = 
(q — x)^/2a -I- h, with the random hypersurface <l?o(q), as the latter is both x and a independent. 
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weighted sum of single-particle velocities ([75|), with the weights given by eq. (|6^, with js from eq. ( |73| ) and Sd from 
eq. 

At early times, the hypersurface "l?ci is a narrow, distorted paraboloid whose apex starts in q = x and gradually 
moves away from it. At these early times, the curvature of this random hypersurface is large, so that there is a 
one-to-one correspondence between Lagrangian and Eulerian space. As time passes, our classical velocity potential 
becomes shallower, so that it becomes increasingly frequent to find points x such that the first-touching condition 
of E/i and $ci occurs at multiple saddle points qs^ . The Lagrangian to Eulerian correspondence has now become 
many-to-one, with the degree of the mapping being inversely proportional to the dimensionality of the structure 
formed. Two simultaneous points of touching, qsi(x, a) and qs2(x, a), signal the formation of a sheet, or pancake, in 
the corresponding Eulerian location x, three indicate the presence of a filament, and four that of a knot. As is well 
known from the adhesion model, filaments occur at the intersection of pancakes, while nodes at the intersection of 
filaments. 

There is only a caveat in our geometrical construction, owing to the approximate form of eq. ([7l|). Because of 
this, eq. ( |65[ ) is only approximately satisfied, namely Vq$ci calculated from eq. ( [7l| ) contains spurious terms which are 
of second-order in the displacement vector. As a consequence, our accuracy in the determination of the saddle points 
is limited to first order. Nevertheless, our Lagrangian approximation is likely to be quite accurate not only at early 
times, when the evolution is linear, but also at later times, as far as the skeleton of the IGM structure is concerned. 
In fact, the dominant contribution to our saddle-point solution in denser regions, such as pancakes, filaments and 
knots, is expected to come from the deepest potential wells [see, e.g., the related interpretation of the path-integral 
solution by Zel'dovich et al. (1985, 1987)], i.e. precisely from those Lagrangian sites where the particles experience 
the smallest displacement from their initial positions (e.g. Matarrese et al. 1992). 

As for the standard adhesion model, the above geometrical construction cannot describe the internal structure 
of dense regions, like pancakes, filaments and knots, as it is only designed to find the sites where these structure 
form and to follow the subsequent evolution through a continuous merging process (e.g. Shandarin & Zel'dovich 
1989). From this point of view, it is clear that a numerical algorithm based on the direct integration of the solution 
in eq. (|6l]), would be far more suitable to follow the baryon distribution on smaller scales. This kind of numerical 
integration can be performed by extending the technique used by Weinberg & Gunn ( 1990a, b) and Melott, Shandarin 
Weinberg (1994). For the free adhesion approximation, Weinberg & Gunn notice that the exact solution [our eq. 



(B3)] appears as a Gaussian convolution integral, which is easy to handle numerically. For our stochastic adhesion 
model, this type of calculation is obviously made more complex by the presence of the random force term, which does 
not generally allow to get an exact explicit expression for the action. This problem, however, is enormously simplified 
if one adopt the approximate expression of eq. (^l|), for the classical action. As mentioned previously, the last step in 
this scheme would consist in the numerical integration of eq. (|6^) to obtain the actual particle positions. Although 
the computational time required to perform the integration in eq. ( |6l] ) plus that of eq. (|6^) implies only a modest 
speed up, compared to a particle-mesh N-body code; one should keep in mind that our stochastic adhesion model 
actually aims at describing a physical situation where the exact treatment of the two-fluid evolution would require a 
full hydrodynamical code. 



7 STATISTICS OF THE DENSITY FIELD 

The stochastic adhesion model provides a direct insight on the properties of the velocity field away from the linear 
regime. A vast literature is devoted to the study of the statistical properties of this field, in cases when the external 
potential is a Gaussian random field which is either time- independent (e.g. Zel'dovich et al. 1987), or with white- 
noise time-correlation properties (e.g. Bouchaud et al. 1995; E et al. 1997). It is generally found that the system 
develops intermittency at late times. Based on this property, Jones (1999) argues that the velocity field is normally 
distributed even during the nonlinear stage. He also suggests that this may provide a dynamical motivation for a 
Lognormal PDF of the density field in the baryonic component. 

More complex is to analyze the statistical behavior of the velocity and density field in our model, owing to the 
non-trivial time-correlation properties of the noise, as well as to the fact that we are interested in the behavior of 
the system in the intermediate regime, when most of the matter still resides in pancakes and filaments, rather than 
having relaxed in the most massive clumps, corresponding to the nodes of the cellular structure. We will only give 
here a preliminary introduction to this problem, focusing on the density field, rather than the peculiar velocity. 

A convenient form of the comoving local matter density can be obtained by integrating the continuity equation 
along the particle trajectories. One gets (e.g. Matarrese et al. 1992) 
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1 + (5(x(q, a), a) — exp 



dr 6'(x(q, r),r) 



(76) 



where S(x, a) = Vx • u(x, a) and the time integration is along the particle trajectory. The validity of this integral 
expression for the density is restricted to the single-stream regime. We are then allowed to apply it in our case, 
provided the viscosity coefficient v takes a finite value. For each Eulerian position x at time a there is, in fact, a 
unique Lagrangian element q implicitly defined through eq. (|6^. 

A relevant statistical information on the density field could then be obtained by studying its (disconnected) 
moments, ((1 + SY). We have 



exp 



-ip-l) [ dr e(x(q, 
Jo 



(77) 



where the subscripts 'E' and 'L' denote the statistical ensemble over which the average is performed: the one on the 
LHS involves the Eulerian PDF, whereas that on the RHS is most conveniently obtained with the Lagrangian PDF. 
To account for this difference we multiplied the RHS by the Jacobian ||9x/9q|| = (1 + 5)~^ (e.g. Kofman et al. 1994). 
Note that the Lagrangian-to-Eulerian mapping is always one-to-one, as the velocity field solving the forced Burgers 
equation is non-singular everywhere for any finite value of v; thus the actual fluid elements, which move according to 
eq. (^), do not experience orbit-crossing during their evolution. 

The main technical difficulty in dealing with eq. ([t^), or any equivalent expression for the density field, resides 
in the inversion of the Lagrangian map q ^ x(q, a) (e.g. Vergassola et al. 1994). The same problem, of course, makes 
it generally difficult to deal analytically with its statistical properties (such as the disconnected moments above). 
An alternative technique to reconstruct the matter density field from the velocity, solving the Burgers equation, is 
discussed by Gurbatov (1996; see also Gurbatov, Malakhov & Saichev 1991), who propose a modification of the 
continuity equation to simplify this inversion problem. Another interesting approach to the same problem is adopted 
by Vergassola et al. (1994), who implement a Fast Legendre Transform algorithm. Ifere we will sketch the first steps 
of a different method, which aims at relating the mass density in x at time a to the properties of the classical velocity 
fleld at the saddle-points qs(x, a). 



7.1 The velocity divergence 

The first step to obtain the density field is to write a suitable expression for the velocity divergence 9. Using eq. (pl[), 
we can write the velocity divergence in the general form 

Next, we can replace for u and U the expressions of eq. (^l|) and (|9|) (which would be exact in the free adhesion 
case), which gives 

e(x,a) = ^-i-[[(uc,(x,q,a)- [[uc,(x,q,a)]])']] , (79) 
where we introduced the symbol 

fd^g A(x,q,a) e-*'->("''''"'/^" 
A(x,q,a) = ^ ^ \ , (80) 

for any scalar, vector or tensor A. To get the first term on the RHS of eq. (j79|) we used the approximate form for the 
classical action deriving from eq. ([7l|). 

To better understand the meaning of the expression (^9|) for the velocity divergence, we can consider the small 
V limit, using the steepest descent approximation. We get 

^(x, a) = ^u)s(qs,T) 6lci(x(qs,a),a) - ^(5u^(x,a) , (81) 
where 6'ci = VxUd = Vx'l'ci, and 

(5u^(x,a) = ^«;s(qs,a)Uci(qa,a) - ^ ^ (q^, a)'!i^s' (q^' , a)uci(qs, a)uci(qs' , a) , (82) 

s s s' 

from which it becomes evident its role as mean square velocity dispersion caused by the convergence of several classical 
streams in x. Note that, by definition, 5v? > 0. In order to obtain the expression of eq. (^l|), we expanded the classical 
velocity to first order around the saddle point. 



© 0000 RAS, MNRAS 000, 000-000 



22 S. Matarrese and R. Mohayaee 



Mci(x,q,a) = <i(x,qs,a) 



and wrote the classical velocity divergence in the form 



6'ci(x(q, a),a) = 



9^-0(q,a) 
dq^dqJ 



dqidqj 



(83) 



(84) 



which follows from taking the Eulerian divergence of eq. (^). Note that, while in the free adhesion case, eq. ( |8l| ) 
is the direct result of the saddle-point approximation, in our stochastic adhesion model, that expression has been 
derived only to first order in the displacement vector. Nonetheless, we believe that the validity of eq. (^ij) is much 
more general than its approximate derivation might suggest. 



7.2 Evolution of the density field 

We are now in a position to discuss the evolution of the density field in the various stages of the structure formation 
process. 



7.2.1 Linear stage 

At early times, when the system is well approximated by Eulerian first-order perturbation theory, drO « —aV^ip ^ 
1, and the density fluctuation is linearly related to that of the initial gravitational potential, which we assumed to be 
a uniform Gaussian process. Density fiuctuations at this time are also Gaussian distributed. 



7.2.2 Weakly nonlinear stage 

Later on, the system enters a weakly nonlinear regime, when density fiuctuations are no longer small, but caustic 
formation is still a sporadic event. In this laminar stage the evolution of the system is generally well described by 
Lagrangian first-order perturbation theory (see Section 4.1). Our solution (^) yields 



1 + (5(x(q,a),a) = exp l^-y dr 6'ci(x(qs, r), r)J , (85) 

as, for each x, there is a unique steepest-descent solution, corresponding to the particle, with initial coordinate qs, 
which reaches that point at time a. This expression is of course equivalent to the more conventional one of Section 
4.1.2 [e.g. Appendix A in (Matarrese et al. 1992)]. At this early nonlinear stage, the density distribution is no longer 
Gaussian. Obtaining the PDF of the mass density at this stage is a classical problem which has been solved long 
ago by Doroshkevich (1970), and analyzed in more detail by several authors (e.g. Kofman et al. 1994; Bernardeau & 
Kofman 1994), for the coUisionless case. As anticipated in Section 4.1.2, the presence of non-zero acceleration in the 
motion of our coUisional fiuid elements does not imply relevant modifications compared with the DM case, so, for 
instance, the functional form of the PDF is unchanged. 



7.2.3 Mildly nonlinear stage 

Let us finally come to the most interesting issue: the analysis of the mildly nonlinear stage, when dense, spatially 
extended structures have grown in the system. It is this process which is accompanied by the occurrence of strong 
phase coherence both in the velocity and density fields. According to our stochastic adhesion model, the formation 
of a shock singularity in x at time a is accompanied by the occurrence of multiple saddle-point solutions qs(x, a). 
As a consequence of these multiple solutions, the velocity dispersion term of eq. ( p2[ ) will suddenly become non-zero. 
Because of the 'diverging' (as i/ —> 0) factor 1/2;^, however, any non-zero Su^ would signal a 'singularity' (for v ^ 0) 
in the local density. One has, therefore, the following picture: as u —> the mass tends to get concentrated in the 
cellular structure, which becomes a thin cobweb of delta-like sheets, filaments and nodes. On the other hand, as soon 
as the Lagrangian-to-Eulerian mapping becomes many-to-one our expression (|76| ) for the mass density will lose its 
validity, and a more clever technique should be used to account for the actual mass converging into these nonlinear 
structures. A possible hint in this direction might be given by the calculation of the velocity field on a pancake, 
made by Kofman (1989), within the free adhesion model. A more detailed analysis of this problem will be presented 
elsewhere. 
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8 CONCLUSIONS 

In this paper we have proposed a new dynamical model for the formation of mildly nonlinear structure in the IGM, 
based on a 'stochastic adhesion' approximation to the Navier-Stokes equation, which governs the evolution of the 
collisional baryon component. The main idea of our model is that a random external force can be used to approximate 
the composition of Hubble drag, local gravitational force and pressure gradients. Such a force, consistently calculated 
from first-order perturbation theory, has two main effects. First, it produces a slight distortion of the large-scale 
filamentary cosmic web that characterizes the IGM distribution, compared with the DM one. Second, it adds a small- 
scale noise to the otherwise deterministic evolution of the fluid, thereby imprinting characteristic ripples on top of 
the cellular network. This small-scale roughness, which is the memory of acoustic oscillations in the baryon fluid, can 
be used as a probe of the IGM equation of state and thermal history after hydrogen reionization. 

This paper was mostly devoted to the presentation of the stochastic adhesion model. A few applications were 
shown, only for the simplified inviscid {u — 0) case, where the trajectories followed by patches of the collisional fluid 
are obtained through a flrst-order Lagrangian approximation scheme, somehow similar to that adopted by Gnedin 
& Hui (1996) and Hui, Gnedin & Zhang (1997). Already at this level, however, the precise IGM thermal history is 
found to affect the resulting clustering properties of the gas, through the detailed structure of denser regions. The 
second part of our paper was devoted to the solution of the forced Burgers equation, which governs the evolution of 
the IGM peculiar velocity field in our model. Among our main results is the saddle-point solution, eq. (^^, and the 
Lagrangian approximation of Section 6: both are completely new and might have much wider applications than the 
present context. We also sketched a geometrical representation of our saddle-point solution, analogous to that used 
for the standard adhesion model, which allows to draw the skeleton of the IGM distribution. 

We did not address here the important issue of whether pancakes or filaments provide the dominant structure 
in the large-scale matter distribution. According to the cosmic web description of Bond et al. (1996) the nonlinear 
evolution of the DM brings a filamentary network into relief. Bond & Wadsley (1997) applied similar ideas to the IGM 
distribution probed by the Lya forest. Quite recently, Colombi, Pogosyan & Souradeep (2000) adopted topological 
descriptors to conclude that the structure of the Universe at the percolating threshold is predominantly filamentary, 
which confirms the results of a percolation analysis of N-body simulations by Sathyaprakash, Sahni & Shandarin 
(1996). In both the free and stochastic adhesion models, this aspect of the emerging structures on large scales is 
expected to be strictly related to the stage of dynamical evolution of the system, as well as to the value of the 
viscosity coefficient u and to the resolution scale at which the structure is analy 

The formation of a cellular structure on large scales, appearing as a cobweb of interconnected sheets, filaments 
and nodes, both for the DM and IGM components, is accompanied by the growth of specific phase coherence in the 
velocity and density field. We have shown how this phenomenon, which has been dubbed 'intermediate intermittency' 
(Zel'dovich et al. 1985, 1987), is most naturally described in terms of the forced Burgers equation, which underlies our 
model. A wavelet analysis of the Lya forest, such as that recently initiated by Meiksin (2000) and Theuns & Zaroubi 
(2000), would be the ideal tool for a statistical analysis of this phenomenon, even on those small scales where the 
effects of the thermal gas pressure are expected to imprint characteristic features on the baryon distribution, which 
should tell us about the equation of state and detailed thermal history of the IGM. 

In concluding this paper let us emphasize that the stochastic adhesion model presented here might represent 
the ideal tool to study the statistical properties of the low-column density Ly-alpha forest, because of its ability to 
accurately describe the IGM distribution from the largest scales, poorly probed by present-day hydro-simulations, 

down to those regions which have undergone mildly non-linear evolution (baryon density contrast up to 5 10), 

which cannot be properly represented by smoothed versions of the Zel'dovich approximation. 
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APPENDIX A: FIRST-ORDER LAGRANGIAN PERTURBATION THEORY FOR THE 
BARYONS 

We will briefly derive here an expression for the baryon trajectories within first-order Lagrangian perturbation theory. 
Lagrangian approximation schemes have been successfully applied to the weakly nonlinear evolution of the DM fluid 
(e.g. Buchert 1989, 1992; Moutarde et al. 1991; Bouchet et al. 1992; Catelan 1995). Adler & Buchert (1999) have 
recently applied the Lagrangian perturbation technique to the case of a collisional self-gravitating fluid, corresponding 
to /b = 1 in our equation (Q) . The case discussed here is instead that of a collisional fluid moving in the gravitational 
field caused by a DM component, corresponding to /b = in our equation (Q). 

The particle trajectories of both components can be represented, before shell-crossing, in the general form 



x(q,a) = q + S(q,a) 

Mass conservation (starting from a uniform distribution) implies 



(Al) 



1 + (5(x(q, a), a) = 



-1 = ' 



V, ■ S + 0(S'") 



Aet{dxa/dqi3) 

where J is the Jacobian determinant of the transformation from Lagrangian to Eulerian coordinates. 
Taking the divergence of the baryon Euler equation we obtain 



Vx-Sb = - — 
la 



Vx- St 



5dm 
a 



1 



(7- l)afc; 



-Vi(l + 5b) 



7-1 



Substituting for the overdensities from eq. (A2) we have 
3 



Vx • Sb + — Vx • Sb 

2a 



Jb 



2a2 Jdm 2(7 - I)a2fc2 ' " V Jb 
Next, we expand everything to first order in the displacement vector S and obtain 



(A2) 



(A3) 



(A4) 
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^•^-^^^ ■ " ■ = ^^V^Vq • Sb . (A5) 

Note that at a given time, the same Eulerian position x is generally reached by the two components (DM and baryons) 
starting from different Lagrangian positions qdm and qb- To lowest order in the displacement vector, however, we 
are allowed to ignore this difference and set simply q = qDM = qb- 

We now assume that the flow is irrotational Vx x S = 0, which, to lowest order implies S = Vq"!/, and reduces 
( |A5| ) to the scalar differential equation 

which can be solved in Fourier space. Using the well-known DM solution \['dm = —aipo, we obtain 

which can be solved in exactly the same manner as we previously solved (^) to obtain the baryon overdensity. The 
solution can be expressed in the form 'i'h = —aiph, where ?/;b(k, a) = Wh{k, a)ipo{'k) with Wh the IGM linear filter 
function obtained in the main text. 



APPENDIX B: PATH-INTEGRAL SOLUTION FOR THE FREE ADHESION MODEL 

In the case of zero potential 77, eq. (^) reduces to the Burgers equation. By dropping the potential rj we obtain the 
free particle trajectory 

x(q,a) = q + auo(q) (Bl) 
and the classical action reduces to 

5.(x,q,a) = = ^ = (^i-^ . (B2) 



2a 

" U 

We therefore have 

W(x,a) = F(a) J d\ e-(Sci(x.q,a)+*o(q))/2. ^ j d\ e-^^-*"''')/^'' , (B3) 

where the pre-factor is evaluated by performing the Gaussian integration over ^ in eq. (^sj). 

The above expression (B3) is equivalent to the standard expression for the adhesion approximation (e.g. Shandarin 
& Zel'dovich 1989). 

In the steepest descent approximation, valid in the limit of vanishing viscosity, the particle trajectories are given 
by the solution of 

Vq ($0 + ScO = , (B4) 
that is 

Vq$o(q) + "^—^ = , (B5) 



which coincides with eq. (Bl), for our set of initial conditions. The results of this Appendix also show that the 
quadratic approximation is exact for free particles. 
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